library(tidyverse) # data manipulation
library(ggpubr) # producing data exploratory plots
library(modelsummary) # descriptive data
library(glmmTMB) # running generalised mixed models
library(DHARMa) # model diagnostics
library(performance) # model diagnostics
library(ggeffects) # partial effect plots
library(car) # running Anova on model
library(emmeans) # post-hoc analysis m1 <- read_csv("import_data/1_month_size_data_2022_2023.csv") |>
mutate(across(1:15,factor)) |>
mutate(STANDARD_LENGTH =LENGTH,
.keep = "unused") |>
select(!(NOTES)) |>
select(1:15,"STANDARD_LENGTH","MASS")
m2 <- read_csv("import_data/2_month_size_data_2022_2023.csv") |>
mutate(across(1:15,factor)) |>
mutate(STANDARD_LENGTH =LENGTH,
.keep = "unused") |>
select(!(NOTES)) |>
select(1:15,"STANDARD_LENGTH","MASS")
m2.5 <- read_csv("import_data/2-5_month_size_data_2022_2023.csv") |>
mutate(across(1:15,factor)) |>
mutate(STANDARD_LENGTH =LENGTH,
.keep = "unused") |>
select(!(NOTES)) |>
select(1:15,"STANDARD_LENGTH","MASS")
adult <- read_csv("import_data/adult_size_2022_2023.csv") |>
mutate(across(1:3,factor),
MALE = FISH_ID,
FEMALE = FISH_ID,
POPULATION = str_sub(FISH_ID, 2,4),
POPULATION = case_when(POPULATION == "ARL" ~ "Arlington Reef",
POPULATION == "SUD" ~ "Sudbury Reef",
POPULATION == "VLA" ~ "Vlassof cay",
POPULATION == "PRE" ~ "Pretty patches",
TRUE ~ POPULATION)) |>
left_join(select(m1, c("MALE","TEMPERATURE")),
by="MALE") |>
left_join(select(m1, c("FEMALE","TEMPERATURE")),
by="FEMALE") |>
distinct() |>
mutate(TEMPERATURE = coalesce(TEMPERATURE.x, TEMPERATURE.y)) |>
drop_na(TEMPERATURE) |>
select(-c("TEMPERATURE.x","TEMPERATURE.y"))m1_df <- m1 |>
left_join(select(adult, c("MALE", "SL", "MASS")),
by ="MALE") |>
mutate(SL_MALE =SL,
MASS_MALE =MASS.y,
.keep = "unused") |>
left_join(select(adult, c("FEMALE", "SL", "MASS")),
by ="FEMALE") |>
mutate(SL_FEMALE =SL,
MASS_FEMALE =MASS,
.keep ="unused") |>
mutate(SL_MIDPOINT = (SL_MALE+SL_FEMALE)/2,
MASS_MIDPOINT = (MASS_MALE+MASS_FEMALE)/2)
m2_df <- m2 |>
left_join(select(adult, c("MALE", "SL", "MASS")),
by ="MALE") |>
mutate(SL_MALE =SL,
MASS_MALE =MASS.y,
.keep = "unused") |>
left_join(select(adult, c("FEMALE", "SL", "MASS")),
by ="FEMALE") |>
mutate(SL_FEMALE =SL,
MASS_FEMALE =MASS,
.keep ="unused") |>
mutate(SL_MIDPOINT = (SL_MALE+SL_FEMALE)/2,
MASS_MIDPOINT = (MASS_MALE+MASS_FEMALE)/2)
m2.5_df <- m2.5 |>
left_join(select(adult, c("MALE", "SL", "MASS")),
by ="MALE") |>
mutate(SL_MALE =SL,
MASS_MALE =MASS.y,
.keep = "unused") |>
left_join(select(adult, c("FEMALE", "SL", "MASS")),
by ="FEMALE") |>
mutate(SL_FEMALE =SL,
MASS_FEMALE =MASS,
.keep ="unused") |>
mutate(SL_MIDPOINT = (SL_MALE+SL_FEMALE)/2,
MASS_MIDPOINT = (MASS_MALE+MASS_FEMALE)/2) |>
drop_na(SL_MALE)plot1 <- ggplot(m1_df, aes(x=SL_MALE, y=STANDARD_LENGTH, color=TEMPERATURE)) +
geom_point(alpha=0.05) +
stat_smooth(method = "lm") +
theme_classic()
plot2 <- ggplot(m1_df, aes(x=SL_FEMALE, y=STANDARD_LENGTH, color=TEMPERATURE)) +
geom_point(alpha=0.05) +
stat_smooth(method = "lm") +
theme_classic()
plot3 <- ggplot(m1_df, aes(x=SL_MIDPOINT, y=STANDARD_LENGTH, color=TEMPERATURE)) +
geom_point(alpha=0.05) +
stat_smooth(method = "lm") +
theme_classic()
ggarrange(plot1, plot2, plot3,
nrow =1,
ncol =3,
common.legend = TRUE)plot1 <- ggplot(m1_df, aes(x=MASS_MALE, y=MASS.x, color=TEMPERATURE)) +
geom_point(alpha=0.05) +
stat_smooth(method = "lm") +
ylim(0,0.15) +
theme_classic()
plot2 <- ggplot(m1_df, aes(x=MASS_FEMALE, y=MASS.x, color=TEMPERATURE)) +
geom_point(alpha=0.05) +
stat_smooth(method = "lm") +
ylim(0,0.15) +
theme_classic()
plot3 <- ggplot(m1_df, aes(x=MASS_MIDPOINT, y=MASS.x, color=TEMPERATURE)) +
geom_point(alpha=0.05) +
stat_smooth(method = "lm") +
ylim(0,0.15) +
theme_classic()
ggarrange(plot1, plot2, plot3,
nrow =1,
ncol =3,
common.legend = TRUE)| POPULATION | 27 | 28.5 | 30 |
|---|---|---|---|
| Arlington Reef | 8 | 8 | 8 |
| Pretty patches | 4 | 6 | 6 |
| Sudbury Reef | 4 | 4 | 2 |
| Vlassof cay | 6 | 2 | 6 |
| POPULATION | 27 | 28.5 | 30 |
|---|---|---|---|
| Arlington Reef | 231 | 105 | 202 |
| Pretty Patches | 116 | 82 | 142 |
| Sudbury Reef | 117 | 55 | 47 |
| Vlassof Cay | 114 | 60 | 120 |
| POPULATION | 27 | 28.5 | 30 |
|---|---|---|---|
| Arlington Reef | 198 | 108 | 152 |
| Pretty Patches | 113 | 83 | 134 |
| Sudbury Reef | 113 | 60 | 56 |
| Vlassof Cay | 77 | 58 | 113 |
| POPULATION | 27 | 28.5 | 30 |
|---|---|---|---|
| Arlington Reef | 239 | 103 | 211 |
| Pretty Patches | 100 | 83 | 133 |
| Sudbury Reef | 111 | 53 | 50 |
| Vlassof Cay | 102 | 60 | 106 |
datasummary(Factor(TEMPERATURE) ~ SL * (NUnique + mean + median + min + max + sd + Histogram),
data = drop_na(adult, SL),
fmt = "%.2f")| TEMPERATURE | NUnique | mean | median | min | max | sd | Histogram |
|---|---|---|---|---|---|---|---|
| 27 | 21 | 97.14 | 100.60 | 84.76 | 104.42 | 7.00 | ▃▂▁▁▂▁▆▇ |
| 28.5 | 20 | 91.72 | 93.83 | 72.47 | 100.00 | 7.57 | ▁▁▂▃▅▇▃ |
| 30 | 22 | 92.52 | 93.22 | 80.34 | 101.22 | 6.44 | ▂▅▂▃▂▃▇▂▇▅ |
datasummary(Factor(TEMPERATURE) ~ STANDARD_LENGTH * (NUnique + mean + median + min + max + sd + Histogram),
data = drop_na(m1_df, STANDARD_LENGTH),
fmt = "%.2f")| TEMPERATURE | NUnique | mean | median | min | max | sd | Histogram |
|---|---|---|---|---|---|---|---|
| 27 | 402 | 12.94 | 12.99 | 8.25 | 18.05 | 1.86 | ▁▂▃▆▇▇▄▃▁ |
| 28.5 | 220 | 13.26 | 13.23 | 8.74 | 17.87 | 1.45 | ▂▄▇▅▆▁ |
| 30 | 344 | 13.51 | 13.48 | 7.73 | 17.76 | 1.63 | ▁▃▅▇▆▄▂ |
datasummary(Factor(TEMPERATURE) ~ STANDARD_LENGTH * (NUnique + mean + median + min + max + sd + Histogram),
data = drop_na(m2_df, STANDARD_LENGTH),
fmt = "%.2f")| TEMPERATURE | NUnique | mean | median | min | max | sd | Histogram |
|---|---|---|---|---|---|---|---|
| 27 | 368 | 20.36 | 20.81 | 10.73 | 27.53 | 2.41 | ▁▁▂▆▇▃ |
| 28.5 | 264 | 20.27 | 20.47 | 12.87 | 26.83 | 2.40 | ▁▁▄▅▇▆▂▁ |
| 30 | 351 | 20.08 | 20.28 | 11.69 | 26.14 | 2.58 | ▁▁▁▁▄▇▇▄▂ |
datasummary(Factor(TEMPERATURE) ~ STANDARD_LENGTH * (NUnique + mean + median + min + max + sd + Histogram),
data = drop_na(m2.5_df, STANDARD_LENGTH),
fmt = "%.2f")| TEMPERATURE | NUnique | mean | median | min | max | sd | Histogram |
|---|---|---|---|---|---|---|---|
| 27 | 439 | 24.15 | 24.56 | 11.92 | 31.70 | 3.02 | ▁▁▂▄▇▅▂ |
| 28.5 | 267 | 23.65 | 23.98 | 14.69 | 30.15 | 2.98 | ▁▁▁▃▅▆▇▆▃▁ |
| 30 | 382 | 23.58 | 23.94 | 13.74 | 31.04 | 2.61 | ▁▁▂▄▇▇▃▁ |
datasummary(Factor(TEMPERATURE) ~ MASS * (NUnique + mean + median + min + max + sd + Histogram),
data = drop_na(adult, MASS),
fmt = "%.2f")| TEMPERATURE | NUnique | mean | median | min | max | sd | Histogram |
|---|---|---|---|---|---|---|---|
| 27 | 21 | 46.93 | 50.63 | 29.85 | 59.93 | 10.30 | ▇▁▃▄▄▆▄ |
| 28.5 | 20 | 38.81 | 42.36 | 16.28 | 53.26 | 10.26 | ▁▁▄▃▄▇▆▁ |
| 30 | 22 | 39.94 | 39.62 | 23.91 | 57.31 | 9.43 | ▅▇▂▇▇▇▂▇▅▂ |
datasummary(Factor(TEMPERATURE) ~ MASS.x * (NUnique + mean + median + min + max + sd + Histogram),
data = drop_na(m1_df, MASS.x),
fmt = "%.2f")| TEMPERATURE | NUnique | mean | median | min | max | sd | Histogram |
|---|---|---|---|---|---|---|---|
| 27 | 437 | 0.06 | 0.06 | 0.01 | 0.22 | 0.03 | ▃▇▇▄▂▁ |
| 28.5 | 248 | 0.07 | 0.06 | 0.02 | 0.19 | 0.02 | ▁▅▇▄▃ |
| 30 | 389 | 0.07 | 0.07 | 0.01 | 0.15 | 0.02 | ▄▅▇▆▅▃▂▁ |
datasummary(Factor(TEMPERATURE) ~ MASS.x * (NUnique + mean + median + min + max + sd + Histogram),
data = drop_na(m2_df, MASS.x),
fmt = "%.2f")| TEMPERATURE | NUnique | mean | median | min | max | sd | Histogram |
|---|---|---|---|---|---|---|---|
| 27 | 460 | 0.27 | 0.27 | 0.03 | 0.68 | 0.10 | ▁▂▄▇▆▃▁ |
| 28.5 | 297 | 0.28 | 0.27 | 0.06 | 0.64 | 0.10 | ▁▃▇▇▇▄▁▁▁ |
| 30 | 425 | 0.28 | 0.27 | 0.04 | 0.62 | 0.11 | ▂▂▅▇▇▄▂▁▁ |
datasummary(Factor(TEMPERATURE) ~ MASS.x * (NUnique + mean + median + min + max + sd + Histogram),
data = drop_na(m2.5_df, MASS.x),
fmt = "%.2f")| TEMPERATURE | NUnique | mean | median | min | max | sd | Histogram |
|---|---|---|---|---|---|---|---|
| 27 | 529 | 0.46 | 0.46 | 0.04 | 1.00 | 0.17 | ▁▂▄▇▇▆▄▁▁ |
| 28.5 | 293 | 0.45 | 0.44 | 0.08 | 1.00 | 0.18 | ▂▄▆▇▆▅▃▂▁ |
| 30 | 477 | 0.45 | 0.45 | 0.08 | 0.96 | 0.15 | ▁▂▄▆▇▄▂▁ |
modelNULL <- glmmTMB(STANDARD_LENGTH ~ 1,
family=gaussian(),
data =m1_df)
model1 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER),
family=gaussian(),
data =m1_df)
model2 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER) + (1|LEVEL),
family=gaussian(),
data = m1_df)
model3 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER) + (1|CLUTCH_ORDER),
family=gaussian(),
data = m1_df)
model4 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER) + (1|REGION),
family=gaussian(),
data = m1_df)
model5 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER) + (1|REGION) + (1|POPULATION),
family=gaussian(),
data = m1_df)
model6 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER) + (1|CLUTCH_ORDER) + (1|REGION) + (1|POPULATION),
family=gaussian(),
data = m1_df)
AIC(modelNULL, model1, model2, model3, model4, model5, model6) modelNULL <- glmmTMB(STANDARD_LENGTH ~ 1,
family=gaussian(),
data =m2_df)
model1 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER),
family=gaussian(),
data =m2_df)
model2 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER) + (1|LEVEL),
family=gaussian(),
data = m2_df)
model3 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER) + (1|CLUTCH_ORDER),
family=gaussian(),
data = m2_df)
model4 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER) + (1|REGION),
family=gaussian(),
data = m2_df)
model5 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER) + (1|REGION) + (1|POPULATION),
family=gaussian(),
data = m2_df)
model6 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER) + (1|CLUTCH_ORDER) + (1|REGION) + (1|POPULATION),
family=gaussian(),
data = m2_df)
AIC(modelNULL, model1, model2, model3, model4, model5, model6) modelNULL <- glmmTMB(STANDARD_LENGTH ~ 1,
family=gaussian(),
data =m2.5_df)
model1 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER),
family=gaussian(),
data =m2.5_df)
model2 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER) + (1|LEVEL),
family=gaussian(),
data = m2.5_df)
model3 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER) + (1|CLUTCH_ORDER),
family=gaussian(),
data = m2.5_df)
model4 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER) + (1|REGION),
family=gaussian(),
data = m2.5_df)
model5 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER) + (1|REGION) + (1|POPULATION),
family=gaussian(),
data = m2.5_df)
model6 <- glmmTMB(STANDARD_LENGTH ~ (1|CLUTCH_NUMBER) + (1|CLUTCH_ORDER) + (1|REGION) + (1|POPULATION),
family=gaussian(),
data = m2.5_df)
AIC(modelNULL, model1, model2, model3, model4, model5, model6) For standard length measurements at differrent time periods, including 1, 2, and 2.5 months the best model is the most simply model (i.e., model1), where the only random factor that is present is CLUTCH_NUMBER.
Now that we have figured out which random factors will be included within out generalized linear mixed effects model we can start to explore different hypothesese by adding in our fixed factors - covariates.
Fixed factors that will be included will be those that are essential to answering the initial research question based on heiritability of traits between offspring and parental fish - labelled as MALE and FEMALE in the dataframe as well as their combined score MIDPOINT, if applicable. TEMPERATURE is also essential to answering the main research question that looks to see if heritability changes at different temperatures.
Our main research hypothesis will be modelled using the formula below”
An alternative research hypothesis will will test will include an interaction with PARENTAL_DAYS_IN_TEMPERATURE to see if heritability was affect by how long adults spent at experimental temperatures. This model may look something like:
Lets start fitting models:
## Warning in AIC.default(modelNULL, model1a, model1b, k = 12): models are not all
## fitted to the same number of observations
## Warning in BIC.default(modelNULL, model1a, model1b): models are not all fitted
## to the same number of observations
Model1a was selected as the best model and will be used going forward.
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.876 0.92 0.98 0.568 0.712 0.62 0.916 0.7 0.296 0.888 0.196 0.596 0.612 0.812 0.536 0.22 0.884 0.62 0.52 0.924 ...
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.021579, p-value = 0.5377
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.013, p-value = 0.832
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 13, observations = 1388, p-value = 0.5431
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.004996154 0.015962883
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.009365994
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.021579, p-value = 0.5377
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.013, p-value = 0.832
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 13, observations = 1388, p-value = 0.5431
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.004996154 0.015962883
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.009365994
## `check_outliers()` does not yet support models of class `glmmTMB`.
## NOTE: Results may be misleading due to involvement in interactions
## Family: gaussian ( identity )
## Formula: STANDARD_LENGTH ~ SL_MALE * TEMPERATURE + (1 | CLUTCH_NUMBER)
## Data: m1_df
##
## AIC BIC logLik deviance df.resid
## 4522.3 4564.2 -2253.2 4506.3 1380
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## CLUTCH_NUMBER (Intercept) 1.204 1.097
## Residual 1.339 1.157
## Number of obs: 1388, groups: CLUTCH_NUMBER, 50
##
## Dispersion estimate for gaussian family (sigma^2): 1.34
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.705421 3.835079 0.445 0.65654
## SL_MALE 0.115953 0.039683 2.922 0.00348 **
## TEMPERATURE28.5 1.173998 6.946791 0.169 0.86580
## TEMPERATURE30 9.270902 5.175407 1.791 0.07324 .
## SL_MALE:TEMPERATURE28.5 -0.003166 0.074194 -0.043 0.96596
## SL_MALE:TEMPERATURE30 -0.088563 0.055022 -1.610 0.10749
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2.5 % 97.5 % Estimate
## (Intercept) -5.81119565 9.22203719 1.705420772
## SL_MALE 0.03817567 0.19373130 0.115953484
## TEMPERATURE28.5 -12.44146230 14.78945897 1.173998335
## TEMPERATURE30 -0.87270857 19.41451247 9.270901952
## SL_MALE:TEMPERATURE28.5 -0.14858388 0.14225099 -0.003166446
## SL_MALE:TEMPERATURE30 -0.19640438 0.01927813 -0.088563126
## Std.Dev.(Intercept)|CLUTCH_NUMBER 0.89522033 1.34545606 1.097487865
## # R2 for Mixed Models
##
## Conditional R2: 0.537
## Marginal R2: 0.120
m1.sl <- emmeans(model1a, ~ SL_MALE*TEMPERATURE,
at =list(SL_MALE=seq(from =min(m1_df$SL_MALE), to =max(m1_df$SL_MALE), by=.25)))
m1.sl.df <- as.data.frame(m1.sl)
m1.sl.obs <- drop_na(m1_df, SL_MALE, STANDARD_LENGTH) |>
mutate(Pred =predict(model1a, re.form=NA, type ='response'),
Resid =residuals(model1a, type ='response'),
Fit =Pred+Resid)
m1.sl.obs.summarize <- m1.sl.obs |>
group_by(CLUTCH_NUMBER, TEMPERATURE) |>
summarise(mean.sl =mean(STANDARD_LENGTH, na.rm=TRUE),
mean.sl.male =mean(SL_MALE, na.rm = TRUE),
sd.sl =sd(STANDARD_LENGTH, na.rm =TRUE),
n.sl = n()) |>
mutate(se.sl = sd.sl / sqrt(n.sl),
lower.ci.sl =mean.sl - qt(1-(0.05/2), n.sl -1) * se.sl,
upper.ci.sl =mean.sl + qt(1-(0.05/2), n.sl -1) * se.sl) |>
ungroup()## `summarise()` has grouped output by 'CLUTCH_NUMBER'. You can override using the
## `.groups` argument.
ggplot(data = m1.sl.df, aes(x=SL_MALE, y=emmean)) +
stat_smooth(aes(color=TEMPERATURE),
method = "lm") +
geom_pointrange(data = m1.sl.obs.summarize, aes(x =mean.sl.male,
y =mean.sl,
ymin =lower.ci.sl,
ymax =upper.ci.sl,
color = TEMPERATURE)) +
scale_color_manual(values = c("#69d7d8","#ff9c56", "#903146")) +
scale_fill_manual(values =c("#69d7d8","#ff9c56", "#903146")) +
facet_wrap(~TEMPERATURE)+
xlab("PARENTAL MALE STANDARD LENGTH (mm)") +
ylab("OFFSPRING STANDARD LENGTH (mm)") +
ggtitle("Offspring-male relationship") +
theme_classic()## `geom_smooth()` using formula = 'y ~ x'
## Warning in AIC.default(modelNULL, model1a, model1b, k = 12): models are not all
## fitted to the same number of observations
## Warning in BIC.default(modelNULL, model1a, model1b): models are not all fitted
## to the same number of observations
The null model appears better than the models that we used. Let’s explore the data bit more and see if we can find a reason for this. Let’s start by looking at a basic histogram of our data.
There appears to be a left skew within our data. Let’s see if this can be better modelled with a Gamma distribution. If not we can try to incorporate transformations to our response variable. The model validations below also could use some improving.
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0 0.124 0.028 0.068 0.4 0.384 0.396 0.492 0.476 0.328 0.504 0.776 0.664 0.612 0.184 0.548 0.644 0.408 0.384 0.588 ...
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.069304, p-value = 1.066e-05
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.99972, p-value = 1
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 24, observations = 1264, p-value =
## 0.0001642
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.01220249 0.02812061
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.01898734
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.069304, p-value = 1.066e-05
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.99972, p-value = 1
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 24, observations = 1264, p-value =
## 0.0001642
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.01220249 0.02812061
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.01898734
## `check_outliers()` does not yet support models of class `glmmTMB`.
model1a.gammaid <- glmmTMB(STANDARD_LENGTH ~ SL_MALE*TEMPERATURE + (1|CLUTCH_NUMBER),
family=Gamma(link="identity"),
data=m2_df)## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in AIC.default(modelNULL, model1a, model1a.gammalog, model1a.gammaid, :
## models are not all fitted to the same number of observations
## Warning in BIC.default(modelNULL, model1a, model1a.gammalog, model1a.gammaid):
## models are not all fitted to the same number of observations
Conclusion: Using a gamma distribution does not help the model. Next we can try transforming the data. A sqrt transformation will be used on the response variable to help with the left skewness.
model1a.sqrt <- glmmTMB(sqrt(STANDARD_LENGTH) ~ SL_MALE*TEMPERATURE + (1|CLUTCH_NUMBER),
family=gaussian(link="identity"),
data=m2_df)
model1a.loglink <- glmmTMB(sqrt(STANDARD_LENGTH) ~ SL_MALE*TEMPERATURE + (1|CLUTCH_NUMBER),
family=gaussian(link="log"),
data=m2_df) ## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in AIC.default(modelNULL, model1a, model1a.sqrt, model1a.loglink, :
## models are not all fitted to the same number of observations
## Warning in BIC.default(modelNULL, model1a, model1a.sqrt, model1a.loglink, :
## models are not all fitted to the same number of observations
Conclusion: The log transformation greatly improves model performance. The sqrt-transformed model will be used moving forward. Lets move on to model validations
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0 0.132 0.028 0.068 0.42 0.432 0.416 0.508 0.496 0.34 0.528 0.756 0.668 0.62 0.196 0.56 0.648 0.42 0.396 0.616 ...
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.098532, p-value = 4.387e-11
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.99991, p-value = 0.984
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 30, observations = 1264, p-value =
## 2.488e-07
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.01606931 0.03370971
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.02373418
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.098532, p-value = 4.387e-11
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.99991, p-value = 0.984
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 30, observations = 1264, p-value =
## 2.488e-07
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.01606931 0.03370971
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.02373418
## `check_outliers()` does not yet support models of class `glmmTMB`.
## Model has sqrt-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the sqrt-scale.
## NOTE: Results may be misleading due to involvement in interactions
## Model has sqrt-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the sqrt-scale.
## Family: gaussian ( identity )
## Formula:
## sqrt(STANDARD_LENGTH) ~ SL_MALE * TEMPERATURE + (1 | CLUTCH_NUMBER)
## Data: m2_df
##
## AIC BIC logLik deviance df.resid
## 336.3 377.5 -160.2 320.3 1256
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## CLUTCH_NUMBER (Intercept) 0.008566 0.09255
## Residual 0.071513 0.26742
## Number of obs: 1264, groups: CLUTCH_NUMBER, 47
##
## Dispersion estimate for gaussian family (sigma^2): 0.0715
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.6342522 0.3940090 11.762 <2e-16 ***
## SL_MALE -0.0013400 0.0040633 -0.330 0.742
## TEMPERATURE28.5 -0.1533775 0.6713293 -0.228 0.819
## TEMPERATURE30 -0.0046237 0.5333702 -0.009 0.993
## SL_MALE:TEMPERATURE28.5 0.0014952 0.0071495 0.209 0.834
## SL_MALE:TEMPERATURE30 -0.0003783 0.0056418 -0.067 0.947
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2.5 % 97.5 % Estimate
## (Intercept) 3.862008774 5.406495540 4.6342521569
## SL_MALE -0.009303963 0.006623978 -0.0013399925
## TEMPERATURE28.5 -1.469158754 1.162403822 -0.1533774658
## TEMPERATURE30 -1.050010054 1.040762723 -0.0046236656
## SL_MALE:TEMPERATURE28.5 -0.012517595 0.015507955 0.0014951799
## SL_MALE:TEMPERATURE30 -0.011436096 0.010679509 -0.0003782931
## Std.Dev.(Intercept)|CLUTCH_NUMBER 0.070908531 0.120808073 0.0925544328
## # R2 for Mixed Models
##
## Conditional R2: 0.110
## Marginal R2: 0.003
m2.sl <- emmeans(model1a.sqrt, ~ SL_MALE*TEMPERATURE,
at =list(SL_MALE=seq(from =min(m2_df$SL_MALE), to =max(m2_df$SL_MALE), by=.25)))
m2.sl.df <- as.data.frame(m2.sl)
m2.sl.obs <- drop_na(m2_df, SL_MALE, STANDARD_LENGTH) |>
mutate(Pred =predict(model1a.sqrt, re.form=NA, type ='response'),
Resid =residuals(model1a.sqrt, type ='response'),
Fit =Pred+Resid)
m2.sl.obs.summarize <- m2.sl.obs |>
group_by(CLUTCH_NUMBER, TEMPERATURE) |>
summarise(mean.sl =mean(STANDARD_LENGTH, na.rm=TRUE),
mean.sl.male =mean(SL_MALE, na.rm = TRUE),
sd.sl =sd(STANDARD_LENGTH, na.rm =TRUE),
n.sl = n()) |>
mutate(se.sl = sd.sl / sqrt(n.sl),
lower.ci.sl =mean.sl - qt(1-(0.05/2), n.sl -1) * se.sl,
upper.ci.sl =mean.sl + qt(1-(0.05/2), n.sl -1) * se.sl) |>
ungroup()## `summarise()` has grouped output by 'CLUTCH_NUMBER'. You can override using the
## `.groups` argument.
ggplot(data = m2.sl.df, aes(x=SL_MALE, y=emmean)) +
stat_smooth(aes(color=TEMPERATURE),
method = "lm",
formula = y^2 ~ x) +
geom_pointrange(data = m2.sl.obs.summarize, aes(x =mean.sl.male,
y =mean.sl,
ymin =lower.ci.sl,
ymax =upper.ci.sl,
color = TEMPERATURE)) +
scale_color_manual(values = c("#69d7d8","#ff9c56", "#903146")) +
scale_fill_manual(values =c("#69d7d8","#ff9c56", "#903146")) +
facet_wrap(~TEMPERATURE)+
xlab("PARENTAL MALE STANDARD LENGTH (mm)") +
ylab("OFFSPRING STANDARD LENGTH (mm)") +
ggtitle("Offspring-male relationship") +
theme_classic()Once again the NULL model seems to outperform our hypothesis testing models. Let’s follow the steps that we conducted for 2-month data and appy a log transformation to our dataset to see if it improved the model.
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.944 0.876 0.86 0.564 0.652 0.788 0.86 0.42 0.896 0.404 0.94 0.272 0.788 0.824 0.928 0.712 0.892 0.972 0.848 0.336 ...
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.078695, p-value = 1.194e-07
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98858, p-value = 0.792
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 32, observations = 1343, p-value =
## 9.318e-08
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.01635359 0.03347159
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.02382725
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.078695, p-value = 1.194e-07
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98858, p-value = 0.792
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 32, observations = 1343, p-value =
## 9.318e-08
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.01635359 0.03347159
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.02382725
## `check_outliers()` does not yet support models of class `glmmTMB`.
## Model has log-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the log-scale.
## NOTE: Results may be misleading due to involvement in interactions
## Model has log-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the log-scale.
## Family: gaussian ( identity )
## Formula:
## log(STANDARD_LENGTH) ~ SL_MALE * TEMPERATURE + (1 | CLUTCH_NUMBER)
## Data: m2.5_df
##
## AIC BIC logLik deviance df.resid
## -1775.2 -1733.6 895.6 -1791.2 1335
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## CLUTCH_NUMBER (Intercept) 0.002184 0.04673
## Residual 0.014515 0.12048
## Number of obs: 1343, groups: CLUTCH_NUMBER, 52
##
## Dispersion estimate for gaussian family (sigma^2): 0.0145
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.001805 0.179768 16.698 <2e-16 ***
## SL_MALE 0.001841 0.001861 0.989 0.322
## TEMPERATURE28.5 0.091645 0.322762 0.284 0.776
## TEMPERATURE30 0.187153 0.238310 0.785 0.432
## SL_MALE:TEMPERATURE28.5 -0.001078 0.003446 -0.313 0.754
## SL_MALE:TEMPERATURE30 -0.002231 0.002536 -0.880 0.379
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2.5 % 97.5 % Estimate
## (Intercept) 2.649466952 3.354143381 3.001805167
## SL_MALE -0.001806051 0.005487995 0.001840972
## TEMPERATURE28.5 -0.540956737 0.724247427 0.091645345
## TEMPERATURE30 -0.279926559 0.654232769 0.187153105
## SL_MALE:TEMPERATURE28.5 -0.007832090 0.005676656 -0.001077717
## SL_MALE:TEMPERATURE30 -0.007201632 0.002739709 -0.002230961
## Std.Dev.(Intercept)|CLUTCH_NUMBER 0.036225309 0.060280991 0.046730050
## # R2 for Mixed Models
##
## Conditional R2: 0.141
## Marginal R2: 0.011
m1.sl <- emmeans(model1a.log, ~ SL_MALE*TEMPERATURE,
at =list(SL_MALE=seq(from =min(m2.5_df$SL_MALE), to =max(m2.5_df$SL_MALE), by=.25)))
m1.sl.df <- as.data.frame(m1.sl)
m1.sl.obs <- drop_na(m2.5_df, SL_MALE, STANDARD_LENGTH) |>
mutate(Pred =predict(model1a.log, re.form=NA, type ='response'),
Resid =residuals(model1a.log, type ='response'),
Fit =Pred+Resid)
m1.sl.obs.summarize <- m1.sl.obs |>
group_by(CLUTCH_NUMBER, TEMPERATURE) |>
summarise(mean.sl =mean(Fit, na.rm=TRUE),
mean.sl.male =mean(SL_MALE, na.rm = TRUE),
sd.sl =sd(Fit, na.rm =TRUE),
n.sl = n()) |>
mutate(se.sl = sd.sl / sqrt(n.sl),
lower.ci.sl =mean.sl - qt(1-(0.05/2), n.sl -1) * se.sl,
upper.ci.sl =mean.sl + qt(1-(0.05/2), n.sl -1) * se.sl) |>
ungroup()## `summarise()` has grouped output by 'CLUTCH_NUMBER'. You can override using the
## `.groups` argument.
ggplot(data = m1.sl.df, aes(x=SL_MALE, y=emmean)) +
stat_smooth(aes(color=TEMPERATURE),
method = "lm",
formula = y ~ x) +
geom_pointrange(data = m1.sl.obs.summarize, aes(x =mean.sl.male,
y =mean.sl,
ymin =lower.ci.sl,
ymax =upper.ci.sl,
color = TEMPERATURE)) +
scale_color_manual(values = c("#69d7d8","#ff9c56", "#903146")) +
scale_fill_manual(values =c("#69d7d8","#ff9c56", "#903146")) +
facet_wrap(~TEMPERATURE)+
xlab("PARENTAL MALE STANDARD LENGTH (mm)") +
ylab("OFFSPRING STANDARD LENGTH (mm)") +
ggtitle("Offspring-male relationship") +
theme_classic()## Warning in AIC.default(modelNULL, model1a, model1b, k = 12): models are not all
## fitted to the same number of observations
## Warning in BIC.default(modelNULL, model1a, model1b): models are not all fitted
## to the same number of observations
Model1a was selected as the best model and will be used going forward.
## qu = 0.25, log(sigma) = -2.465176 : outer Newton did not converge fully.
## qu = 0.25, log(sigma) = -2.373479 : outer Newton did not converge fully.
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.884 0.944 0.996 0.688 0.8 0.732 0.964 0.74 0.364 0.952 0.304 0.724 0.78 0.86 0.6 0.348 0.912 0.768 0.544 0.936 ...
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.018029, p-value = 0.7799
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0184, p-value = 0.808
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 11, observations = 1331, p-value = 0.8766
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.004132601 0.014739191
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.008264463
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.018029, p-value = 0.7799
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0184, p-value = 0.808
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 11, observations = 1331, p-value = 0.8766
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.004132601 0.014739191
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.008264463
## `check_outliers()` does not yet support models of class `glmmTMB`.
## NOTE: Results may be misleading due to involvement in interactions
## Family: gaussian ( identity )
## Formula:
## STANDARD_LENGTH ~ SL_FEMALE * TEMPERATURE + (1 | CLUTCH_NUMBER)
## Data: m1_df
##
## AIC BIC logLik deviance df.resid
## 4340.7 4382.2 -2162.3 4324.7 1323
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## CLUTCH_NUMBER (Intercept) 1.099 1.048
## Residual 1.346 1.160
## Number of obs: 1331, groups: CLUTCH_NUMBER, 48
##
## Dispersion estimate for gaussian family (sigma^2): 1.35
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.53924 2.90633 1.906 0.0567 .
## SL_FEMALE 0.07590 0.03017 2.516 0.0119 *
## TEMPERATURE28.5 2.68132 4.55262 0.589 0.5559
## TEMPERATURE30 5.32954 4.88350 1.091 0.2751
## SL_FEMALE:TEMPERATURE28.5 -0.02056 0.04863 -0.423 0.6725
## SL_FEMALE:TEMPERATURE30 -0.04828 0.05141 -0.939 0.3476
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2.5 % 97.5 % Estimate
## (Intercept) -0.15707023 11.2355419 5.53923585
## SL_FEMALE 0.01677235 0.1350219 0.07589712
## TEMPERATURE28.5 -6.24165530 11.6042999 2.68132230
## TEMPERATURE30 -4.24194847 14.9010207 5.32953612
## SL_FEMALE:TEMPERATURE28.5 -0.11587486 0.0747581 -0.02055838
## SL_FEMALE:TEMPERATURE30 -0.14904158 0.0524739 -0.04828384
## Std.Dev.(Intercept)|CLUTCH_NUMBER 0.85086513 1.2919121 1.04844787
## # R2 for Mixed Models
##
## Conditional R2: 0.505
## Marginal R2: 0.101
m1_df <-
m1_df |>
drop_na(SL_FEMALE)
m2.sl <- emmeans(model1a.sqrt, ~ SL_MALE*TEMPERATURE,
at =list(SL_MALE=seq(from =min(m2_df$SL_MALE), to =max(m2_df$SL_MALE), by=.25)))
m1.sl <- emmeans(model1a, ~ SL_FEMALE*TEMPERATURE,
at =list(SL_FEMALE=seq(from =min(m1_df$SL_FEMALE), to =max(m1_df$SL_FEMALE), by=.25)))
m1.sl.df <- as.data.frame(m1.sl)
m1.sl.obs <- drop_na(m1_df, SL_FEMALE, STANDARD_LENGTH) |>
mutate(Pred =predict(model1a, re.form=NA, type ='response'),
Resid =residuals(model1a, type ='response'),
Fit =Pred+Resid)
m1.sl.obs.summarize <- m1.sl.obs |>
group_by(CLUTCH_NUMBER, TEMPERATURE) |>
summarise(mean.sl =mean(Fit, na.rm=TRUE),
mean.sl.female =mean(SL_FEMALE, na.rm = TRUE),
sd.sl =sd(Fit, na.rm =TRUE),
n.sl = n()) |>
mutate(se.sl = sd.sl / sqrt(n.sl),
lower.ci.sl =mean.sl - qt(1-(0.05/2), n.sl -1) * se.sl,
upper.ci.sl =mean.sl + qt(1-(0.05/2), n.sl -1) * se.sl) |>
ungroup()## `summarise()` has grouped output by 'CLUTCH_NUMBER'. You can override using the
## `.groups` argument.
ggplot(data = m1.sl.df, aes(x=SL_FEMALE, y=emmean)) +
stat_smooth(aes(color=TEMPERATURE),
method = "lm") +
geom_pointrange(data = m1.sl.obs.summarize, aes(x =mean.sl.female,
y =mean.sl,
ymin =lower.ci.sl,
ymax =upper.ci.sl,
color = TEMPERATURE)) +
scale_color_manual(values = c("#69d7d8","#ff9c56", "#903146")) +
scale_fill_manual(values =c("#69d7d8","#ff9c56", "#903146")) +
facet_wrap(~TEMPERATURE)+
xlab("PARENTAL FEMALE STANDARD LENGTH (mm)") +
ylab("OFFSPRING STANDARD LENGTH (mm)") +
ggtitle("Offspring-male relationship") +
theme_classic()## `geom_smooth()` using formula = 'y ~ x'
## Warning in AIC.default(modelNULL, model1a, model1b, k = 12): models are not all
## fitted to the same number of observations
## Warning in BIC.default(modelNULL, model1a, model1b): models are not all fitted
## to the same number of observations
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.004 0.092 0.024 0.096 0.404 0.328 0.392 0.512 0.412 0.32 0.52 0.78 0.6 0.572 0.172 0.464 0.604 0.352 0.324 0.512 ...
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.076927, p-value = 1.165e-06
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0001, p-value = 0.984
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 17, observations = 1213, p-value = 0.03297
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.008184788 0.022344553
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.01401484
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.076927, p-value = 1.165e-06
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0001, p-value = 0.984
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 17, observations = 1213, p-value = 0.03297
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.008184788 0.022344553
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.01401484
## `check_outliers()` does not yet support models of class `glmmTMB`.
model1a.sqrt <- glmmTMB(sqrt(STANDARD_LENGTH) ~ SL_FEMALE*TEMPERATURE + (1|CLUTCH_NUMBER),
family=gaussian(link="identity"),
data=m2_df)
model1a.pwr3 <- glmmTMB(STANDARD_LENGTH^3 ~ SL_FEMALE*TEMPERATURE + (1|CLUTCH_NUMBER),
family=gaussian(link="identity"),
data=m2_df)
model1a.log <- glmmTMB(log(STANDARD_LENGTH) ~ SL_FEMALE*TEMPERATURE + (1|CLUTCH_NUMBER),
family=gaussian(link="identity"),
data=m2_df)## Warning in AIC.default(modelNULL, model1a, model1a.sqrt, model1a.pwr3,
## model1a.log, : models are not all fitted to the same number of observations
## Warning in BIC.default(modelNULL, model1a, model1a.sqrt, model1a.pwr3,
## model1a.log): models are not all fitted to the same number of observations
Conclusion: The log transformation greatly improves model performance. The sqrt-transformed model will be used moving forward. Lets move on to model validations
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0 0.116 0.024 0.104 0.448 0.364 0.416 0.536 0.444 0.34 0.548 0.776 0.62 0.608 0.188 0.48 0.612 0.388 0.368 0.56 ...
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.10529, p-value = 4.171e-12
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0004, p-value = 0.976
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 25, observations = 1213, p-value =
## 2.499e-05
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.01338107 0.03027497
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.02061006
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.10529, p-value = 4.171e-12
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0004, p-value = 0.976
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 25, observations = 1213, p-value =
## 2.499e-05
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.01338107 0.03027497
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.02061006
## `check_outliers()` does not yet support models of class `glmmTMB`.
## Model has sqrt-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the sqrt-scale.
## NOTE: Results may be misleading due to involvement in interactions
## Model has sqrt-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the sqrt-scale.
## Family: gaussian ( identity )
## Formula:
## sqrt(STANDARD_LENGTH) ~ SL_FEMALE * TEMPERATURE + (1 | CLUTCH_NUMBER)
## Data: m2_df
##
## AIC BIC logLik deviance df.resid
## 305.3 346.1 -144.6 289.3 1205
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## CLUTCH_NUMBER (Intercept) 0.007992 0.0894
## Residual 0.070568 0.2656
## Number of obs: 1213, groups: CLUTCH_NUMBER, 45
##
## Dispersion estimate for gaussian family (sigma^2): 0.0706
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.5186835 0.2995276 15.086 <2e-16 ***
## SL_FEMALE -0.0002317 0.0030872 -0.075 0.940
## TEMPERATURE28.5 -0.2971534 0.4491625 -0.662 0.508
## TEMPERATURE30 0.5367662 0.5186774 1.035 0.301
## SL_FEMALE:TEMPERATURE28.5 0.0032236 0.0047758 0.675 0.500
## SL_FEMALE:TEMPERATURE30 -0.0059220 0.0054242 -1.092 0.275
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2.5 % 97.5 % Estimate
## (Intercept) 3.931620211 5.105746782 4.5186834966
## SL_FEMALE -0.006282586 0.005819150 -0.0002317182
## TEMPERATURE28.5 -1.177495697 0.583188921 -0.2971533882
## TEMPERATURE30 -0.479822825 1.553355148 0.5367661612
## SL_FEMALE:TEMPERATURE28.5 -0.006136814 0.012583936 0.0032235609
## SL_FEMALE:TEMPERATURE30 -0.016553152 0.004709128 -0.0059220122
## Std.Dev.(Intercept)|CLUTCH_NUMBER 0.067828174 0.117834224 0.0894006727
## # R2 for Mixed Models
##
## Conditional R2: 0.110
## Marginal R2: 0.009
m2_df <-
m2_df |>
drop_na(SL_FEMALE)
m2.sl <- emmeans(model1a.sqrt, ~ SL_FEMALE*TEMPERATURE,
at =list(SL_FEMALE=seq(from =min(m2_df$SL_FEMALE), to =max(m2_df$SL_FEMALE), by=.25)))
m2.sl.df <- as.data.frame(m2.sl)
m2.sl.obs <- drop_na(m2_df, SL_FEMALE, STANDARD_LENGTH) |>
mutate(Pred =predict(model1a.sqrt, re.form=NA, type ='response'),
Resid =residuals(model1a.sqrt, type ='response'),
Fit =Pred+Resid)
m2.sl.obs.summarize <- m2.sl.obs |>
group_by(CLUTCH_NUMBER, TEMPERATURE) |>
summarise(mean.sl =mean(Fit, na.rm=TRUE),
mean.sl.male =mean(SL_FEMALE, na.rm = TRUE),
sd.sl =sd(Fit, na.rm =TRUE),
n.sl = n()) |>
mutate(se.sl = sd.sl / sqrt(n.sl),
lower.ci.sl =mean.sl - qt(1-(0.05/2), n.sl -1) * se.sl,
upper.ci.sl =mean.sl + qt(1-(0.05/2), n.sl -1) * se.sl) |>
ungroup()## `summarise()` has grouped output by 'CLUTCH_NUMBER'. You can override using the
## `.groups` argument.
ggplot(data = m2.sl.df, aes(x=SL_FEMALE, y=emmean)) +
stat_smooth(aes(color=TEMPERATURE),
method = "lm",
formula = y ~ x) +
geom_pointrange(data = m2.sl.obs.summarize, aes(x =mean.sl.male,
y =mean.sl,
ymin =lower.ci.sl,
ymax =upper.ci.sl,
color = TEMPERATURE)) +
scale_color_manual(values = c("#69d7d8","#ff9c56", "#903146")) +
scale_fill_manual(values =c("#69d7d8","#ff9c56", "#903146")) +
facet_wrap(~TEMPERATURE)+
xlab("PARENTAL FEMALE STANDARD LENGTH (mm)") +
ylab("OFFSPRING STANDARD LENGTH (mm)") +
ggtitle("Offspring-male relationship") +
theme_classic()## Warning in AIC.default(modelNULL, model1a, model1b, k = 12): models are not all
## fitted to the same number of observations
## Warning in BIC.default(modelNULL, model1a, model1b): models are not all fitted
## to the same number of observations
Once again the NULL model seems to outperform our hypothesis testing models. Let’s follow the steps that we conducted for 2-month data and appy a log transformation to our dataset to see if it improved the model.
## Warning in AIC.default(modelNULL, model1a, model1a.log, model1a.sqrt, k = 5):
## models are not all fitted to the same number of observations
## Warning in BIC.default(modelNULL, model1a, model1a.log, model1a.sqrt): models
## are not all fitted to the same number of observations
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.904 0.888 0.812 0.5 0.664 0.764 0.84 0.42 0.84 0.456 0.94 0.352 0.7 0.836 0.928 0.748 0.884 0.98 0.804 0.328 ...
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.07065, p-value = 5.159e-06
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98886, p-value = 0.864
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 25, observations = 1289, p-value =
## 9.858e-05
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.01258969 0.02849825
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.01939488
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.07065, p-value = 5.159e-06
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98886, p-value = 0.864
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 25, observations = 1289, p-value =
## 9.858e-05
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.01258969 0.02849825
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.01939488
## `check_outliers()` does not yet support models of class `glmmTMB`.
## Model has log-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the log-scale.
## NOTE: Results may be misleading due to involvement in interactions
## Model has log-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the log-scale.
## Family: gaussian ( identity )
## Formula:
## log(STANDARD_LENGTH) ~ SL_FEMALE * TEMPERATURE + (1 | CLUTCH_NUMBER)
## Data: m2.5_df
##
## AIC BIC logLik deviance df.resid
## -1722.1 -1680.8 869.0 -1738.1 1281
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## CLUTCH_NUMBER (Intercept) 0.002247 0.04741
## Residual 0.014285 0.11952
## Number of obs: 1289, groups: CLUTCH_NUMBER, 50
##
## Dispersion estimate for gaussian family (sigma^2): 0.0143
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.0169019 0.1436697 20.999 <2e-16 ***
## SL_FEMALE 0.0017149 0.0014920 1.149 0.250
## TEMPERATURE28.5 0.0437181 0.2237656 0.195 0.845
## TEMPERATURE30 0.2064255 0.2355460 0.876 0.381
## SL_FEMALE:TEMPERATURE28.5 -0.0005848 0.0023902 -0.245 0.807
## SL_FEMALE:TEMPERATURE30 -0.0024581 0.0024848 -0.989 0.323
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2.5 % 97.5 % Estimate
## (Intercept) 2.735314415 3.298489462 3.016901939
## SL_FEMALE -0.001209370 0.004639164 0.001714897
## TEMPERATURE28.5 -0.394854344 0.482290509 0.043718083
## TEMPERATURE30 -0.255236128 0.668087148 0.206425510
## SL_FEMALE:TEMPERATURE28.5 -0.005269533 0.004099951 -0.000584791
## SL_FEMALE:TEMPERATURE30 -0.007328298 0.002411998 -0.002458150
## Std.Dev.(Intercept)|CLUTCH_NUMBER 0.036649207 0.061320595 0.047406236
## # R2 for Mixed Models
##
## Conditional R2: 0.149
## Marginal R2: 0.015
m2.5_df <- m2.5_df |>
drop_na(SL_FEMALE)
m2.5.sl <- emmeans(model1a.log, ~ SL_FEMALE*TEMPERATURE,
at =list(SL_FEMALE=seq(from =min(m2.5_df$SL_FEMALE), to =max(m2.5_df$SL_FEMALE), by=.25)))
m2.5.sl.df <- as.data.frame(m2.5.sl)
m2.5.sl.obs <- drop_na(m2.5_df, SL_FEMALE, STANDARD_LENGTH) |>
mutate(Pred =predict(model1a.log, re.form=NA, type ='response'),
Resid =residuals(model1a.log, type ='response'),
Fit =Pred+Resid)
m2.5.sl.obs.summarize <- m2.5.sl.obs |>
group_by(CLUTCH_NUMBER, TEMPERATURE) |>
summarise(mean.sl =mean(Fit, na.rm=TRUE),
mean.sl.male =mean(SL_FEMALE, na.rm = TRUE),
sd.sl =sd(Fit, na.rm =TRUE),
n.sl = n()) |>
mutate(se.sl = sd.sl / sqrt(n.sl),
lower.ci.sl =mean.sl - qt(1-(0.05/2), n.sl -1) * se.sl,
upper.ci.sl =mean.sl + qt(1-(0.05/2), n.sl -1) * se.sl) |>
ungroup()## `summarise()` has grouped output by 'CLUTCH_NUMBER'. You can override using the
## `.groups` argument.
ggplot(data = m2.5.sl.df, aes(x=SL_FEMALE, y=emmean)) +
stat_smooth(aes(color=TEMPERATURE),
method = "lm",
formula = y ~ x) +
geom_pointrange(data = m2.5.sl.obs.summarize, aes(x =mean.sl.male,
y =mean.sl,
ymin =lower.ci.sl,
ymax =upper.ci.sl,
color = TEMPERATURE)) +
scale_color_manual(values = c("#69d7d8","#ff9c56", "#903146")) +
scale_fill_manual(values =c("#69d7d8","#ff9c56", "#903146")) +
facet_wrap(~TEMPERATURE)+
xlab("PARENTAL MALE STANDARD LENGTH (mm)") +
ylab("OFFSPRING STANDARD LENGTH (mm)") +
ggtitle("Offspring-male relationship") +
theme_classic()## Warning in AIC.default(modelNULL, model1a, model1b, k = 12): models are not all
## fitted to the same number of observations
## Warning in BIC.default(modelNULL, model1a, model1b): models are not all fitted
## to the same number of observations
Model1a was selected as the best model and will be used going forward.
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.884 0.94 0.996 0.684 0.8 0.724 0.964 0.732 0.344 0.952 0.28 0.724 0.768 0.856 0.588 0.34 0.912 0.768 0.532 0.936 ...
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.022999, p-value = 0.4821
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0184, p-value = 0.808
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 10, observations = 1331, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.003608554 0.013773413
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.007513148
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.022999, p-value = 0.4821
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0184, p-value = 0.808
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 10, observations = 1331, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.003608554 0.013773413
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.007513148
## `check_outliers()` does not yet support models of class `glmmTMB`.
## NOTE: Results may be misleading due to involvement in interactions
## Family: gaussian ( identity )
## Formula:
## STANDARD_LENGTH ~ SL_MIDPOINT * TEMPERATURE + (1 | CLUTCH_NUMBER)
## Data: m1_df
##
## AIC BIC logLik deviance df.resid
## 4338.3 4379.8 -2161.1 4322.3 1323
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## CLUTCH_NUMBER (Intercept) 1.044 1.022
## Residual 1.346 1.160
## Number of obs: 1331, groups: CLUTCH_NUMBER, 48
##
## Dispersion estimate for gaussian family (sigma^2): 1.35
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.710960 3.212695 1.155 0.24805
## SL_MIDPOINT 0.094769 0.033312 2.845 0.00444 **
## TEMPERATURE28.5 1.737605 5.493089 0.316 0.75176
## TEMPERATURE30 7.100067 4.815692 1.474 0.14038
## SL_MIDPOINT:TEMPERATURE28.5 -0.009484 0.058743 -0.162 0.87174
## SL_MIDPOINT:TEMPERATURE30 -0.066058 0.051056 -1.294 0.19572
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2.5 % 97.5 % Estimate
## (Intercept) -2.58580723 10.00772625 3.710959510
## SL_MIDPOINT 0.02947975 0.16005857 0.094769159
## TEMPERATURE28.5 -9.02865192 12.50386250 1.737605288
## TEMPERATURE30 -2.33851551 16.53864908 7.100066787
## SL_MIDPOINT:TEMPERATURE28.5 -0.12461817 0.10564964 -0.009484266
## SL_MIDPOINT:TEMPERATURE30 -0.16612479 0.03400911 -0.066057843
## Std.Dev.(Intercept)|CLUTCH_NUMBER 0.82884321 1.25961036 1.021772722
## # R2 for Mixed Models
##
## Conditional R2: 0.504
## Marginal R2: 0.120
m1_df <-
m1_df |>
drop_na(SL_MIDPOINT)
m1.sl <- emmeans(model1a, ~ SL_MIDPOINT*TEMPERATURE,
at =list(SL_MIDPOINT=seq(from =min(m1_df$SL_MIDPOINT), to =max(m1_df$SL_MIDPOINT), by=.25)))
m1.sl.df <- as.data.frame(m1.sl)
m1.sl.obs <- drop_na(m1_df, SL_MIDPOINT, STANDARD_LENGTH) |>
mutate(Pred =predict(model1a, re.form=NA, type ='response'),
Resid =residuals(model1a, type ='response'),
Fit =Pred+Resid)
m1.sl.obs.summarize <- m1.sl.obs |>
group_by(CLUTCH_NUMBER, TEMPERATURE) |>
summarise(mean.sl =mean(Fit, na.rm=TRUE),
mean.sl.female =mean(SL_MIDPOINT, na.rm = TRUE),
sd.sl =sd(Fit, na.rm =TRUE),
n.sl = n()) |>
mutate(se.sl = sd.sl / sqrt(n.sl),
lower.ci.sl =mean.sl - qt(1-(0.05/2), n.sl -1) * se.sl,
upper.ci.sl =mean.sl + qt(1-(0.05/2), n.sl -1) * se.sl) |>
ungroup()## `summarise()` has grouped output by 'CLUTCH_NUMBER'. You can override using the
## `.groups` argument.
ggplot(data = m1.sl.df, aes(x=SL_MIDPOINT, y=emmean)) +
stat_smooth(aes(color=TEMPERATURE),
method = "lm") +
geom_pointrange(data = m1.sl.obs.summarize, aes(x =mean.sl.female,
y =mean.sl,
ymin =lower.ci.sl,
ymax =upper.ci.sl,
color = TEMPERATURE)) +
scale_color_manual(values = c("#69d7d8","#ff9c56", "#903146")) +
scale_fill_manual(values =c("#69d7d8","#ff9c56", "#903146")) +
facet_wrap(~TEMPERATURE)+
xlab("PARENTAL MIDPOINT STANDARD LENGTH (mm)") +
ylab("OFFSPRING STANDARD LENGTH (mm)") +
ggtitle("Offspring-male relationship") +
theme_classic()## `geom_smooth()` using formula = 'y ~ x'
## Warning in AIC.default(modelNULL, model1a, model1b, k = 12): models are not all
## fitted to the same number of observations
## Warning in BIC.default(modelNULL, model1a, model1b): models are not all fitted
## to the same number of observations
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.004 0.112 0.024 0.104 0.424 0.348 0.4 0.524 0.436 0.328 0.528 0.792 0.62 0.612 0.184 0.476 0.612 0.368 0.336 0.54 ...
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.075278, p-value = 2.141e-06
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0001, p-value = 0.992
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 17, observations = 1213, p-value = 0.03297
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.008184788 0.022344553
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.01401484
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.075278, p-value = 2.141e-06
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0001, p-value = 0.992
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 17, observations = 1213, p-value = 0.03297
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.008184788 0.022344553
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.01401484
## `check_outliers()` does not yet support models of class `glmmTMB`.
## Warning in AIC.default(modelNULL, model1a, model1a.sqrt, model1a.log, k = 4):
## models are not all fitted to the same number of observations
## Warning in BIC.default(modelNULL, model1a, model1a.sqrt, model1a.log): models
## are not all fitted to the same number of observations
Conclusion: The log transformation greatly improves model performance. The log-transformed model will be used moving forward. Lets move on to model validations
## qu = 0.75, log(sigma) = -2.524445 : outer Newton did not converge fully.
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0 0.128 0.024 0.108 0.452 0.38 0.444 0.556 0.464 0.348 0.556 0.784 0.64 0.62 0.196 0.48 0.636 0.396 0.388 0.56 ...
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.10294, p-value = 1.368e-11
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0004, p-value = 0.984
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 26, observations = 1213, p-value =
## 8.987e-06
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.01404833 0.03124958
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.02143446
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.10294, p-value = 1.368e-11
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0004, p-value = 0.984
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 26, observations = 1213, p-value =
## 8.987e-06
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.01404833 0.03124958
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.02143446
## `check_outliers()` does not yet support models of class `glmmTMB`.
## Model has sqrt-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the sqrt-scale.
## NOTE: Results may be misleading due to involvement in interactions
## Model has sqrt-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the sqrt-scale.
## Family: gaussian ( identity )
## Formula:
## sqrt(STANDARD_LENGTH) ~ SL_MIDPOINT * TEMPERATURE + (1 | CLUTCH_NUMBER)
## Data: m2_df
##
## AIC BIC logLik deviance df.resid
## 306.7 347.5 -145.3 290.7 1205
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## CLUTCH_NUMBER (Intercept) 0.008334 0.09129
## Residual 0.070568 0.26565
## Number of obs: 1213, groups: CLUTCH_NUMBER, 45
##
## Dispersion estimate for gaussian family (sigma^2): 0.0706
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.5680403 0.3463791 13.188 <2e-16 ***
## SL_MIDPOINT -0.0007419 0.0035725 -0.208 0.835
## TEMPERATURE28.5 -0.3097405 0.5622775 -0.551 0.582
## TEMPERATURE30 0.2568173 0.5282549 0.486 0.627
## SL_MIDPOINT:TEMPERATURE28.5 0.0033209 0.0059915 0.554 0.579
## SL_MIDPOINT:TEMPERATURE30 -0.0030421 0.0055681 -0.546 0.585
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2.5 % 97.5 % Estimate
## (Intercept) 3.889149826 5.246930855 4.5680403408
## SL_MIDPOINT -0.007743953 0.006260161 -0.0007418964
## TEMPERATURE28.5 -1.411784137 0.792303094 -0.3097405213
## TEMPERATURE30 -0.778543243 1.292177900 0.2568173288
## SL_MIDPOINT:TEMPERATURE28.5 -0.008422291 0.015064054 0.0033208814
## SL_MIDPOINT:TEMPERATURE30 -0.013955311 0.007871187 -0.0030420618
## Std.Dev.(Intercept)|CLUTCH_NUMBER 0.069457445 0.119981456 0.0912885828
## # R2 for Mixed Models
##
## Conditional R2: 0.110
## Marginal R2: 0.005
m2_df <-
m2_df |>
drop_na(SL_MIDPOINT)
m2.sl <- emmeans(model1a.sqrt, ~ SL_MIDPOINT*TEMPERATURE,
at =list(SL_MIDPOINT=seq(from =min(m2_df$SL_MIDPOINT), to =max(m2_df$SL_MIDPOINT), by=.25)))
m2.sl.df <- as.data.frame(m2.sl)
m2.sl.obs <- drop_na(m2_df, SL_MIDPOINT, STANDARD_LENGTH) |>
mutate(Pred =predict(model1a.sqrt, re.form=NA, type ='response'),
Resid =residuals(model1a.sqrt, type ='response'),
Fit =Pred+Resid)
m2.sl.obs.summarize <- m2.sl.obs |>
group_by(CLUTCH_NUMBER, TEMPERATURE) |>
summarise(mean.sl =mean(Fit, na.rm=TRUE),
mean.sl.male =mean(SL_MIDPOINT, na.rm = TRUE),
sd.sl =sd(Fit, na.rm =TRUE),
n.sl = n()) |>
mutate(se.sl = sd.sl / sqrt(n.sl),
lower.ci.sl =mean.sl - qt(1-(0.05/2), n.sl -1) * se.sl,
upper.ci.sl =mean.sl + qt(1-(0.05/2), n.sl -1) * se.sl) |>
ungroup()## `summarise()` has grouped output by 'CLUTCH_NUMBER'. You can override using the
## `.groups` argument.
ggplot(data = m2.sl.df, aes(x=SL_MIDPOINT, y=emmean)) +
stat_smooth(aes(color=TEMPERATURE),
method = "lm",
formula = y ~ x) +
geom_pointrange(data = m2.sl.obs.summarize, aes(x =mean.sl.male,
y =mean.sl,
ymin =lower.ci.sl,
ymax =upper.ci.sl,
color = TEMPERATURE)) +
scale_color_manual(values = c("#69d7d8","#ff9c56", "#903146")) +
scale_fill_manual(values =c("#69d7d8","#ff9c56", "#903146")) +
facet_wrap(~TEMPERATURE)+
xlab("PARENTAL MIDPOINT STANDARD LENGTH (mm)") +
ylab("OFFSPRING STANDARD LENGTH (mm)") +
ggtitle("Offspring-male relationship") +
theme_classic()## Warning in AIC.default(modelNULL, model1a, model1b, k = 12): models are not all
## fitted to the same number of observations
## Warning in BIC.default(modelNULL, model1a, model1b): models are not all fitted
## to the same number of observations
Once again the NULL model seems to outperform our hypothesis testing models. Let’s follow the steps that we conducted for 2-month data and appy a log transformation to our dataset to see if it improved the model.
## Warning in AIC.default(modelNULL, model1a, model1a.log, model1a.sqrt, k = 5):
## models are not all fitted to the same number of observations
## Warning in BIC.default(modelNULL, model1a, model1a.log, model1a.sqrt): models
## are not all fitted to the same number of observations
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.904 0.888 0.812 0.5 0.664 0.764 0.84 0.42 0.84 0.456 0.94 0.348 0.692 0.832 0.928 0.744 0.88 0.98 0.8 0.328 ...
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.070529, p-value = 5.391e-06
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98904, p-value = 0.864
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 24, observations = 1289, p-value =
## 0.000198
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.01196512 0.02757771
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.01861908
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.070529, p-value = 5.391e-06
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98904, p-value = 0.864
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: simulationOutput
## outliers at both margin(s) = 24, observations = 1289, p-value =
## 0.000198
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.01196512 0.02757771
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.01861908
## `check_outliers()` does not yet support models of class `glmmTMB`.
## Model has log-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the log-scale.
## NOTE: Results may be misleading due to involvement in interactions
## Model has log-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the log-scale.
## Family: gaussian ( identity )
## Formula:
## log(STANDARD_LENGTH) ~ SL_MIDPOINT * TEMPERATURE + (1 | CLUTCH_NUMBER)
## Data: m2.5_df
##
## AIC BIC logLik deviance df.resid
## -1721.8 -1680.5 868.9 -1737.8 1281
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## CLUTCH_NUMBER (Intercept) 0.002266 0.0476
## Residual 0.014284 0.1195
## Number of obs: 1289, groups: CLUTCH_NUMBER, 50
##
## Dispersion estimate for gaussian family (sigma^2): 0.0143
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.0021622 0.1632880 18.386 <2e-16 ***
## SL_MIDPOINT 0.0018651 0.0016939 1.101 0.271
## TEMPERATURE28.5 0.0500368 0.2772671 0.180 0.857
## TEMPERATURE30 0.2036245 0.2394468 0.850 0.395
## SL_MIDPOINT:TEMPERATURE28.5 -0.0006479 0.0029646 -0.219 0.827
## SL_MIDPOINT:TEMPERATURE30 -0.0024310 0.0025427 -0.956 0.339
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 2.5 % 97.5 % Estimate
## (Intercept) 2.682123614 3.322200775 3.0021621947
## SL_MIDPOINT -0.001454866 0.005185117 0.0018651255
## TEMPERATURE28.5 -0.493396679 0.593470289 0.0500368050
## TEMPERATURE30 -0.265682655 0.672931581 0.2036244632
## SL_MIDPOINT:TEMPERATURE28.5 -0.006458364 0.005162652 -0.0006478557
## SL_MIDPOINT:TEMPERATURE30 -0.007414637 0.002552685 -0.0024309759
## Std.Dev.(Intercept)|CLUTCH_NUMBER 0.036826284 0.061530783 0.0476019964
## # R2 for Mixed Models
##
## Conditional R2: 0.149
## Marginal R2: 0.014
m2.5_df <- m2.5_df |>
drop_na(SL_MIDPOINT)
m2.5.sl <- emmeans(model1a.log, ~ SL_MIDPOINT*TEMPERATURE,
at =list(SL_MIDPOINT=seq(from =min(m2.5_df$SL_MIDPOINT), to =max(m2.5_df$SL_MIDPOINT), by=.25)))
m2.5.sl.df <- as.data.frame(m2.5.sl)
m2.5.sl.obs <- drop_na(m2.5_df, SL_MIDPOINT, STANDARD_LENGTH) |>
mutate(Pred =predict(model1a.log, re.form=NA, type ='response'),
Resid =residuals(model1a.log, type ='response'),
Fit =Pred+Resid)
m2.5.sl.obs.summarize <- m2.5.sl.obs |>
group_by(CLUTCH_NUMBER, TEMPERATURE) |>
summarise(mean.sl =mean(Fit, na.rm=TRUE),
mean.sl.male =mean(SL_MIDPOINT, na.rm = TRUE),
sd.sl =sd(Fit, na.rm =TRUE),
n.sl = n()) |>
mutate(se.sl = sd.sl / sqrt(n.sl),
lower.ci.sl =mean.sl - qt(1-(0.05/2), n.sl -1) * se.sl,
upper.ci.sl =mean.sl + qt(1-(0.05/2), n.sl -1) * se.sl) |>
ungroup()## `summarise()` has grouped output by 'CLUTCH_NUMBER'. You can override using the
## `.groups` argument.
ggplot(data = m2.5.sl.df, aes(x=SL_MIDPOINT, y=emmean)) +
stat_smooth(aes(color=TEMPERATURE),
method = "lm",
formula = y ~ x) +
geom_pointrange(data = m2.5.sl.obs.summarize, aes(x =mean.sl.male,
y =mean.sl,
ymin =lower.ci.sl,
ymax =upper.ci.sl,
color = TEMPERATURE)) +
scale_color_manual(values = c("#69d7d8","#ff9c56", "#903146")) +
scale_fill_manual(values =c("#69d7d8","#ff9c56", "#903146")) +
facet_wrap(~TEMPERATURE)+
xlab("PARENTAL MALE STANDARD LENGTH (mm)") +
ylab("OFFSPRING STANDARD LENGTH (mm)") +
ggtitle("Offspring-male relationship") +
theme_classic()